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Abstract 



Detailed ionic models of cardiac cells are difficult for numerical simulations because they consist 
of a large number of equations and contain small parameters. The presence of small parameters, 
however, may be used for asymptotic reduction of the models. Earlier results have shown that 
the asymptotics of cardiac equations are non-standard. Here we apply such a novel asymptotic 
method to an ionic model of human atrial tissue in order to obtain a reduced but accurate 
model for the description of excitation fronts. Numerical simulations of spiral waves in atrial 
tissue show that wave fronts of propagating action potentials break-up and self-terminate. Our 
model, in particular, yields a simple analytical criterion of propagation block, which is similar in 
purpose but completely different in nature to the 'Maxwell rule' in the FitzHugh-Nagumo type 
models. Our new criterion agrees with direct numerical simulations of break-up of re-entrant 
waves. Key words: excitation; conduction; refractoriness; mathematical model 
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1 Introduction 

Refractoriness is a fundamental characteristic of biological excitable media, including cardiac tis- 
sues. The boundary between absolute and relative refractoriness can be defined as the boundary 
between the ability and the inability of the medium to conduct excitation waves 0|. Transient 
conduction block is thought to be a key event in the initiation of re-entrant arrhythmias and in 
the development and the self-perpetuation of atrial and ventricular fibrillation So it 

is important to understand well the immediate causes and conditions of propagation blocks and 
sudden break-ups in such nonstationary regimes. The aim of the present work is to improve 
this understanding via analysis of a mathematical model of human atrial tissue 0| . 

Kohl et al. 0] distinguish two types of single-cell cardiac models: 'membrane potential mod- 
els' and 'ionic current models'. The membrane potential models attempt to represent cellular 
electrical activity by describing, with a minimal number of equations, the spatio-temporal course 
of changes in membrane potential. Their equations are constructed using a dynamical-systems 
arguments to caricature various properties and processes of cardiac function. Examples of this 
type of models start with the mathematical description of heartbeat as a relaxation oscillator 
by van der Pol and van der Mark and continue to play an important role in describing bio- 
physical behaviour with the the most successful one arguably being the FitzHugh-Nagumo 
equations jiol. 



drV = Dd\V + ey {V - - g), 



(1) 



where V and g are dynamical variables corresponding to the action potential and the cardiac 
current gating variables, ey, eg, 7, and /3 are p arameters and D is a diffusion constant. Further 
examples of such models can be found in 3^ 1^. . among others. An attractive feature 



of this approach is that, along with a reasonable description of excitability, threshold, plateau 
and refractoriness, it focusses on generic equations which can often be treated analytically and 
their dynamical properties can be extended and applied to very different physical, chemical or 
biological problems of similar mathematical structure. The main drawback of these models, 
however, is their lack of an explicit correspondence between model components and constituent 
parts of the biological system, e.g. ion channels and transporter proteins. The second type of 
models, the ionic current models, attempt to model action potential (AP) behaviour on the basis 
of ion fluxes in as much detail as possible in order to fit experimental data and predict behaviour 
under previously untested conditions. A major breakthrough in this direction of cell modelling 
was the work of Hodgkin and Huxley representing the first complete quantitative description 
of the giant squid axon. The ionic concept was applied to cardiac cells by Noble 0, 0] 
and there are now ionic models of sinoatrial node pacemaker cells e.g. (l^ . atrial myocytes 
e.g. li^], Purkinje fibres e.g. ventricular myocytes e.g. il, 2^ and cardiac connective tissue 



cells e.g. |2J|. This is only an incomplete list and the collection of available models continues 
to expand. The ionic models have been successfully applied to study various conditions of 
metabolic activity and excitation-contraction coupling, feedback mechanisms, response to drugs, 
etc. For recent reviews of detailed ionic models, their computational aspects and applications 
we refer to the reviews of Kohl et al. and Clayton However, these models are very 

complicated and have to be studied mostly numerically. Their numerical study is aggravated 
by stiffness of the equations, i.e. broad range of characteristic time scales of dynamic variables 
caused by numerous small parameters of the models. 

An attractive compromise is exemplified by the model of Fenton and Karma [2^ . which 
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combines the simplicity of only three differential equations with realistic description of (crudely) 
the AP shape and (rather nicely) the dependence of the AP duration and front propagation speed 
on the diastolic interval, i.e. 'restitution curves'. Unlike the earlier two-component model by 
Aliev and Panfilov [l5| it has a structure similar to that of true ionic models, and its parameters 
have been fitted to mimick properties of selected four detailed ventricular myocyte models. It 
is simpler than later proposed models of the same "intermediate" kind such as |23]. However, 
this deservedly popular model has not been in any way "derived" from any detailed model, so 
it is only reliable within the phenomenology on which it has been validated, i.e. normal or 
premature APs, but not propagation blocks. 

The problem of conditions for propagation has an elegant solution for the FitzHugh-Nagumo 
system Eqs. ^ and its generalizations, within an asymptotic theorv exploiting the difference of 
time scales of different variables, such as eg <S ey in case of Eqs.n|2g]. The answer is formulated 
in terms of the instantaneous values of the slow variables {g in Eqs.^), and claims that excitation 
will propagate if the definite integral of the kinetic term in the right hand side of the equation 
for the fast variable {V in Eqs. H), between the lower and the upper quasi-stationary states, is 
positive [2^. se e eq. 4.5]. This is similar to Maxwell's 'equal areas' rule in the theory of phase 
transitions |30|, see section 9.3]. In case of Eqs. Q , this rule boils down to an inequality for 
the slow variable g: excitation front will propagate if the value of g at it is less than a certain 
g^:. However, FitzHugh-Nagumo type models completely misrepresent the idiosyncratic 'front 
dissipation' scenario by which propagation block happens in the ionic current models 0. The 
reason is that small parameters in such models appear in essentially different ways from the one 
assumed by the standard asymptotic theory 0,0. So, this elegant 'Maxwell rule' solution is 
not applicable to any realistic models. 

We have developed an alternative asymptotic approach based on special mathematical prop- 
erties of the detailed ionic models, not captured by the standard theory |3]- This approach 
demonstrated excellent quantitative accuracy for APs in isolated Noble-1962 model cells js^ ]. 
and correctly, on a qualitative level, described the front dissipation mechanism of break-up of 
re-entrant waves in 6ourtemanche et al. model of human airial tissue, although quantitative 
correspondence with the full model was poor jHH]- In this paper we suggest, for the first time, a 
refined simplified asymptotic model of a cardiac excitation front, which provides numerically ac- 
curate prediction of the front propagation velocity (within 16%) and its profile (within 0.7 mV). 
It also gives an analytical condition for propagation block in a re-entrant wave, expressed as a 
simple inequality involving the slow inactivation gate j of the fast sodium current. The condi- 
tion is in excellent agreement with results of direct numerical simulations of the Courtemanche 
et al. 0] full ionic model of 21 partial differential equations. 

The paper is organised as follows. In ^we introduce the simplified model equations and 
discuss their properties. Analytical solutions are presented in ^for a piecewise linear 'carica- 
ture' version of our simplified model. Accurate numerical results and a two-dimensional test 
are presented in ^ The paper concludes with a discussion of results and questions open for 
future studies in ^ 



Conditions for propagation in atrial tissue 



4 



2 Mathematical formulation of the model equations 
2.1 Asymptotic reduction 



In this section we briefly summarise the asymptotic arguments of |35|] relevant to our present 
purposes. We re-write Courtemanche et al. 0] model in the following one-parameter form: 

OtV = D [dx + JCdx) V - -, 

(fn(V: e) — m) , ^ , - , 
dTm= ^ > , m{V;Q)=M{V)e{V -Vm), 

dTh= hiV;0)=HiV)9iV,-V), 



Btw - 
drd = 



eThiV) 

{w{V) - w) 

{MV) - Og) 
{d{V) - d) 



drU = F{V,...) (2) 

where D is the voltage diffusion constant, e is a small parameter used for the asymptotics, IC 
is the curvature of the propagating front, 9{) is the Heaviside function, is the sum of all 
currents except the fast sodium current Ine, the dynamic variables V, m, h, Ua, Oa and d are 
defined in 0|, U = (j, Oj, . . . , Na^, Kj, . . . is the vector of all other, slower variables, and F is 
the vector of the corresponding right-hand sides. The rationale for this parameterisation is: 

1. The dynamic variables V ^ m, h, Ua, w, Oa, d are 'fast variables', i.e. they change sig- 
nificantly during the upstroke of a typical AP potential, unlike all other variables which 
change only slightly during that period. The relative speed of the dynamical variables 
is estimated by comparing the magnitude of their corresponding 'time scale functions' as 
shown in Fig. ^a). For a system of differential equations dy/dt = F{y) the time scale 
functions are defined as Ti{y) = | ( dFi/ dyi)~^ | , i = 1 . . . and coincide with the functions 
T already present in Eqs. HI 

2. A specific feature of V is that it is fast only because of one of the terms in the right-hand 
side, the large current /Nai whereas other currents are not that large and so do not have 
the large coefficient in front of them. 

3. The fast sodium current /Na is only large during the upstroke of the AP, and not that 
large otherwise as illustrated in Fig. H^d). This is due to the fact that either gate m or 
gate h or both are almost closed outside the upstroke since their quasistationary values 
rn{V) and h{V) are small there as seen in Fig. E^b). Thus in the limit e — > 0, functions 
fn{V) and h{V) have to be considered zero in certain overlapping intervals V G (— oo, Vr 
and V S [Vh,+oo), and Vh < Vm, hence the representations rn{V;0) = M{V)6{V — Vr 
and h{V; 0) = H{V) e{Vh - V). 



mi 
m.) 
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4. The term JCdxV in the first equation represents the effect of the front curvature for waves 
propagating in two or three spatial dimensions. Derivation of this term using asymptotic 
arguments can be found e.g. in |2^. A simple rule-of-thumb way to understand it is this. 
Imagine a circular wave in two spatial dimensions. The diffusion term in the equation for 
V then has the form D (5^ + dy) V = D (9^ + ^dji) V where R is the polar radius. If R 
at the front is large, its instant curvature IC = 1/R changes slowly as the front propagates, 
and can be replaced with a constant for long time intervals. Considering R as a new X 
coordinate, we then get Eqs.HJ 

These aspects, as applied to the fast sodium current, have been shown to be crucial for the 
correct description of the propag ation block a. In particular, it is important that the /.-gate 
is included among the fast variables. The particular importance of h dynamics at the fringe 
of excitability has been noted before, e.g. for the modified Beeler- Renter model (j^ . A more 
detailed discussion of the parameterisation given by Eqs. |21can be found in jHHl- 

A change of variables^ t = e~^T, x = {eD)~^^'^X , k = (eL')^/^/C and subsequently the limit 
e ^ transforms Eqs. 121 into 

dtV = [dl + Kd^) V - C^//Na(l^, h,j), 

dm = {M{V) e{V - V^) - m)/r^(y), 

dth = {H{v) eiVh -V)- h)/Thiv), 

dtUa = {u^{V) - Ua)/TuSy), 

dtw = {w{V)-w)/T^{V), 

dtOa={o-a{V)-Oa)/ToAV), 
dtd={d{V)-d)/Td{V), 

5tU = 0. (3) 

In other words, we consider the fast time scale on which the upstroke of the AP happens, neglect 
the variations of slow variables during this period as well as all transmembrane currents except 
/Na, as they do not make significant contribution during this period and replace rn and h with 
zero when they are small. 

In the resulting Eqs.|31the first three equations for V , m and h form a closed subsystem, the 
following four equations for n^, Oa and can be solved if V{x^ t) is known but do not affect its 
dynamics, and the rest of the equations state that all other variables remain unchanged. Hence 
we concentrate on the first three equations as the system describing propagation of an AP front 
or its failure. The above derivation procedure does not give a precise definition of the functions 
H{y) and M(y), it only requires that these are reasonably close to h[V) and rn{V) for those 
values of V where these functions are not small. Here 'reasonably close' means that replacement 
of h{V) with H{V) e{Vh - V) and m{V) with M {V) e{V - Vm) does not change significantly the 
solutions of interest, i.e. the propagating fronts. We have found that the simplest approximation 
in the form Miy) = 1, H(y) = 1 works well enough. This is demonstrated in tabled where 
various choices of M{V) and HiV) are tested. So, ultimately, we consider the following system 

change of the value of D is equivalent to rescaling of the spatial coordinate, and is not critical to any of 
the questions considered here. In order to operate with dimensional velocity, we assume the value of the diffusion 
coefficient D = 0.03f 25 mm^/ms, as in our earlier publications |35l . l37l . Increase of the diffusion coefficient to, say, 
D = O.lmm^/ms raises the propagation velocity from 0.28mm/ms in Table^to 0.50mm/ms, in full agreement 
e.g. with results of Xie et al. |38| for the same model. 
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where 



dtV = {dl + Kd,) V + I^^iV) jhm\ (4a) 

dth={e{Vh-v)-h)/TH{v), (4b) 

dtm = {e{V - Vm) - m)/Tm{V), (4c) 



In.{V) = gNaiVNa " V), (5a) 
Tk{V) = {ak{V) + (3k{V))~\ k = h,m, (5b) 
ah{V) = 0.135 e-^^+so)/*^-^ e{-V - 40), 
(3f^{V) = (3.56 eO-°^'^^ + 3.1 x 10^ e"-^'^^) e{-V - 40) 
+ e{V + 40) (0.13(1 + e-(V+io.66)/ii.i))-i^ 
_ 0.32(^ + 47.13) 

am[V) - ^ _ g_o,l(y+47.13) ' 

gNa = 7.8, V}va = 67.53, F/, = -66.66, = -32.7. 



All parameters and functions here are defined as in |a] except the new 'gate threshold' parameters 
Vh and Vm which are chosen from the conditions h{Vh) = 1/2 and m^{Vm) = 1/2. As follows 
from the derivation, variable j, the slow inactivation gate of the fast sodium current, acts as a 
parameter of the model. It is the only one of all slow variables included in the vector U that 
affects our fast subsystem. We say that it describes the 'excitability' of the tissue. Notice that 
it is a multiplier of gNa, so a reduced availability of the fast sodium channels, e.g. as under 
tetrodotoxin or arguably in Brugada syndrome can be formally described by a reduced 
value of the parameter j. 

Before proceeding to the analysis of the simplified three- variable model defined by Eqs.|l]we 
wish to demonstrate that it is a good approximation of the full model of 0] both on a qualitative 
and a quantitative level. On the qualitative level, we show that a temporary obstacle leads to 
a dissipation of the front. This is illustrated in Fig. [21 which shows propagation of the AP into 
a region in time and space where the excitability of the tissue is artificially suppressed. The 
sharp wave fronts of the model of Courtemanche et al. [H| as well as of Eqs.EJstop propagating 
and start to spread diffusively once they reach the blocked zone. The propagation does not 
resume after the block is removed. This behaviour is completely different from that of the 
FitzHugh-Nagumo system of Eqs.nin which even though the propagation is blocked for nearly 
the whole duration of the AP, the wave resumes once the block is removed. Tabled illustrates, 
on the quantitative level, the accuracy of Eqs.|l]as an approximation of the full model of 0. 

It is a popular concept going back to classical works [e.g. |^ that the fast activation gate 
m is considered a 'fast variable' and is 'adiabatically eliminated' since most of the time, except 
possibly during a very short transient, it is close to its quasistationary value m ~ m{V). Hence 
the model can be simplified by replacing m with rn(V) and eliminating the equation for m, 

dtV = dlv + l^,9{v-Vm)jh, 

dth= {9{Vh-V)-h)/Th. (6) 
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We have explored this possibihty for the model of Courtemanche et al. [a] in [35|. Eqs.lHlare 
qualitatively correct, i.e. they still show front dissipation on collision with a temporary obstacle, 
but make a large error in the front propagation speed, as demonstrated in table H 



2.2 Travelling waves and reduction to ODE of the three- variable model 

To find out when propagation of excitation is possible in our simplified model and when it will 
be blocked, we study solutions in the form of propagating fronts as well as the conditions of 
existence of such solutions. 

We look for solutions in the form of a front propagating with a constant speed and shape. 
So we use the ansatz F{z) = F{x + ct) for F = V,h, m where 2: = x + ct is a 'travelling wave 
coordinate' and c is the dimensionless wave speed of the front, related to the dimensional speed 
C by c = {e/ny^^C . Then Eqs. 0] reduce to a system of autonomous ordinary differential 
equations, 

V" = {c-k)V'- hT^iV) jhm^ (7a) 
h' = {cTh{V)y\9{Vh-V)-h), (7b) 

m' ={cTra{V))-\e{V-Vra)-7n), (7c) 

where the boundary conditions are given by 

' ^ " " " " " (8a) 

(8b) 
(8c) 

Here Va and K; are the pre- and post-front voltages. 

Eqs. [7| represent a system of fourth order so its general solution depends on four arbitrary 
constants. Together with constants Va, Vi^ and c this makes seven constants to be determined 
from the six boundary conditions in Eqs. |S1 Thus, we should have a one-parameter family 
of solutions, i.e. one of the parameters (V^,V^,c) can be chosen arbitrary from a certain 
range. A natural choice is because the pre-front voltage acts as an initial condition for a 
propagating front in the tissue, and because in our study of the conditions for propagation it is 
most conveniently treated as a parameter rather than as an unknown. 





-00) 




V{+oo) 


= K;, Va<Vh< Vm. < V^, 


h{- 


-00) 


= 1, 


h{+oo) 


= 0, 


m{- 


-00) 


= 0, 




= 1. 



3 Analytical study of the reduced model 
3.1 An exactly solvable caricature model 

The parameter-counting arguments given in the previous section make it plausible that the 
problem defined by Eqs. [3 with boundary conditions of Eqs. [HI has a one-parameter family of 
travelling wave-front solutions. However, the problem is posed in a highly unusual way since 
the asymptotic pre-front and post-front states are not stable isolated equilibria but belong to 
continua of equilibria and thus are only neutrally stable. We are not aware of any general 
theorems that would guarantee existence of solutions of a nonlinear boundary value-eigenvalue 
problem of this kind. For the two-component model of Eqs. El considered in j^a] this worry has 
been alleviated by the fact that there is a 'caricature' model which has the same structure as 
Eqs. ini including the structure and stability of the equilibrium set and which admits an exact 
and exhaustive analytical study [s^]. Fortunately, a similar 'caricature' exists for our present 
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three-variable problem as well. We replace functions It^a.{V), Th{V) and Tm{V) defined in Eqs.ISl 
with constants. The choice of the constants is somewhat arbitrary. We assume that the events 
in the beginning of the interval z G [^,+00), where V is just above V^, are most important 
for the front propagation. So for numerical illustrations we choose the values of constants /Na? 
Th and Tm as the values of the corresponding functions in Eqs. at some fixed value of the 
voltage V. We set the z axis so that V{0) = V^, and then V{(,) = Vm for some ^ > still to be 
determined. We demand that the solutions for the unknowns V, h and m are continuous and 
that V is smooth at the internal boundary points. 

In this formulation, Eqs. I7bl and Ec] decouple from Eq. EH and from each other and solved 
separately. The solutions of these first-order linear ODE with constant coefficients are given 
by Eqs. Ilflb and llOb . respectively. It follows that in the interval V < Vm, Eq. [23 is a linear 
homogeneous ODE with constant coefficients and its solution given at the first row of Eq. IlUb 
satisfies the boundary conditions V{—oo) = Va, V{0) = and V{(,) = Vm provided that 
the internal boundary point is given by Eq. 1121 To solve the linear inhomogeneous Eq. EH 
in the interval V > Vm we note that its inhomogeneous term / = I^s_{V) j hm^ is a sum of 
exponentials 



n=0 ^ ^ 



(9) 



Bn = 1 



Tm. + nTh 
ThTm 



and terms proportional to riTh will appear in the solution due to the expression for i?„. Imposing 
the boundary conditions at the internal point V{S,) = Vm and at infinity V{oo) = V^, we obtain 
the solution in this interval given at the second row of Eq. llOb . Finally, the wave speed c is fixed 
by Eq. Illbl from the requirement that the solution for V{z) is smooth at the internal boundary 
point ^. To summarise, the solution of Eqs. [7|and|Hlis 



Viz) 



Hz) 



3 

n=0 



-z/{cTh) 



2 < 0, 

z>0, 



(10a) 



(10b) 



m{z) 



0, 



z>i, 



(lOc) 



1 _ e(C-z)/(cT™,)^ 

where the pre-front voltage Vq, the post- front voltage V^j and the wave speed c are related by 



Vui = Vm + InhJ [cTh Tm) 6 



2p-«/(cr,.) 



E 

n=0 



Tm + riTh 



= (c - K){ym - Va) - iNaj CThTm e" 



</icTh) 



E 

n=0 



a„ c 



(11a) 
(lib) 
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the distance between points V = and V = Vm is 




) 



(12) 



and An{c, z) and a„(c) are abbreviations for 



An{c,z) 



an{c) ( n^Th - {Tm + nTh)z 
exp 

Tm + nTh V CThTm 



) 



(13a) 



a„(c) 



nj C(C - k) Th Tm + Tm + nTh 



(13b) 



In the hmit r, 
expected. 



m 



this solution tends to the solution of the two-component model of [4^, as 



The accurate expression in Eq. [S^ for the sodium current /Na(^) vanishes for V = V^a 
which, in particular, means that the transmembrane voltage never exceeds V/va- So, replacing 
this function with a constant changes the properties of the system qualitatively. Even bigger 
discrepancies are expected to occur from replacing the Th{V) and Tm{V) by constants because 
these functions vary by an order of magnitude in the range between the pre- and the post- 
front voltage. It is surprising, however, that even this rough approximation produces results 
which, with exception of the post-front voltage, are within several percent from the solution 
of the detailed ionic model i and certainly capture its qualitative features as can be seen in 
Fig. 131 where the constants are chosen at V = Vm, i.e. lNa.{Vm), Th{Vm) and TmiVm)- This 
relatively good agreement is not due to this special choice of parameter values. Indeed, the 
caricature model and its solution Eqs. EH involve the parameters iNa, Th, Tm, n, Va and j. 
The dependence on the curvature k is negligible in comparison to the deviation of the solution 
Eqs. EH of the caricature model from the numerical solution of the three- variable model Eqs. [3 
The dependence on the pre-front voltage Va and the excitability parameter j is discussed in 
section ini21 and represented in Figs. 0] and El The parameters /Na, Th, Tm, on the other hand, 
are somewhat arbitrary but in order to achieve a good agreement with the original system 
given by Eqs. El we choose these values as the values of the corresponding functions in Eqs. 
at various values of V. In Fig. |llthe relationship between the wave speed c and the excitation 
parameter j for several such choices of V is presented. It can be seen that such a variation of 
the values of /Naj Th, Tm does not lead to significant qualitative changes in the solution Eas. IIUI 
of the caricature model. Figs. 01 and 01 also show, for comparison, the numerical solutions of the 
detailed ionic model of 0| and of the full three- variable model of Eqs. I3 which will be described 
in detail in the next section. 

3.2 The condition for propagation 

Equation lllbl defines c as a smooth function of the parameters within a certain domain. The 
boundary of this domain is associated with the propagation failure. Not all parameters, /Na; 
Th, Tm, K, Va and j, entering Eq. lllbl are of equal importance. We consider here k = and 
postpone the investigation of the effects of curvature to the next section. Parameters /Naj Th and 
Tm represent well-defined properties of the tissue, albeit changeable depending on physiological 
conditions. On the other hand, parameters j and Va are not model constants, but 'slowly 
varying' dynamic quantities: j remains approximately constant throughout the front, and Va 
represents the transmembrane voltage ahead of the front, but both can vary widely on large 
scales between different fronts. Hence we need to determine the singular points of the dispersion 
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relation in Eq. Illbl with respect to j and Va ■ 

Similarly to the two-component caricature [s^l, Eg. Illbl is a transcendental equation for c, 
but it is easily solvable for the excitation parameter j: 



{Vrn- Vg) ^ 



e --h JJ (c^ Th Tm + T^ + n Th) • (14) 



The resulting relationship of j and c for a selected value of Va is shown in Fig. 0] This figure 
reveals a bifurcation. For values of j lower than some jmin no travelling wave solutions exist. 
After a bifurcation at j > jmm two solutions with different speeds are possible. Our direct 
numerical simulations of Eqs. as well as studies of the two-component caricature model by 
Hinch suggest that the solutions of the lower branch are unstable. The bifurcation point 
imin can be determined from the condition that j has a minimum with respect to c at this point 
and therefore satisfies 

ocj Va=const 

This produces, with j(c) defined by Eg. 1141 a quintic polynomial equation for c^. 

Activation of the sodium current is possible because Tm Th, permitting transient channel 
opening and current flow through the cell membrane. The ratio Th/Tm is a function of V 
in the full model, and is a constant in Eqs. d The minimal value of this ratio, necessary for 
propagation, is shown on Fig.lSJas a function of various choices of /Nai Tm and j; it is obtained by 
numerical solution of the algebraic equation Eg. Illbl The smallness of Tm/Th allows approximate 
solution of the above mentioned quintic equation for c^. We set 



x;5„c. (16) 



n=0 



Substituting this expansion in Eq. ^jand discarding the small terms of order 0{Tm) gives the 
zeroth-order approximation to the solution as a function of the pre-front voltage Va 

.(0)^ = (K^^^) geWe^+4e (g + 2 + + 49) , (17) 

Q = \u{{Vm-Va)l{Vh-Va)). 



This limit corresponds to the two- variable caricature ^1||. For any given value of the pre- 
front voltage the value of j must be larger than jmin in order for wave fronts to propagate. 
Although lacking sufficient accuracy, the zeroth-order approximation given by Eq. ll7l reproduces 
qualitatively well the conditions for propagation and dissipation of excitation fronts in the model 
of Courtemanche et al. 0. Analogously, discarding the small terms of order O(t^) gives the 
first-order approximation, 

£1 = e-T^^ HiATm - n A), (18) 

O ^ Tm „ 

71=0 



A = 12t|V0(0 + 4), 

A = (Q\e + 4) + e^/\e + 2)V@Ttj (ii ^m-Q^ 



Conditions for propagation in atrial tissue 



11 



This approximation is already very good and changes insignificantly as more terms are consid- 
ered in Eq. ^1 see Fig. ^ 

4 Numerical results 

4.1 Propagating front solutions 

We solved Eqs. I31H1 numerically, using the method described in Appendix ^ The results are 
shown in figures 01 11] and [Tj Figure |21 offers a comparison of the shapes of the solution of Eqs.[7| 
with a snapshot of a travelling wave solution of the full model of Courtemanche et al. 0| . The 
values of the wave speed and the post-front voltage are presented in table H and also show an 
excellent agreement. This confirms our assumptions that the fronts of travelling waves in the 
full model have constant speed and shape and thus satisfy an ODE system, and that j remains 
approximately constant during the front. Figure [7| shows the wave speed c as a function of two 
of the parameters of the problem, the pre- front voltage Va and the excitability parameter j. For 
every value of j and Va from a certain domain, two values of the wave speed c are possible, which 
is similar to the solutions of the caricature model. The smaller values of c are not observed in 
the PDE simulation of the full model. This is a strong indication that they are unstable. 

4.2 The condition for propagation 

In this subsection, we report numerical values for the threshold of excitability jmin below which 
wave fronts are not sustainable and have to dissipate, as predicted by the reduced three- variable 
model of Eqs. EHHl Figure |H1 presents jmin as a function of the pre-front voltage Va- The curve 
jminiVa) represents a boundary in the space of the slow variables {V,j) which separates the 
region of relative refractoriness where excitation fronts are possible, even though possibly slowed 
down, from the region of absolute refractoriness where excitation fronts cannot propagate at all. 
In practice, however, we can reduce the condition of the absolute refractoriness even further. 
This is possible because typical APs have their tails very closely following one path on the 
{V,j) plane. This property is known for cardiac models; e.g. [s^ presents an evidence for the 
Modified Beeler-Reuter model that the dynamics of recovery from an AP do not depend on 
details of how that AP has been initiated. Therefore of the whole curve {V, j^iniV)) only one 
point is important, its intersection with the curve {V{t),j(t)) representing a typical AP tail. 
For the model Courtemanche et al. 0| considered here, we simply state the existence of this 
universal (y{t),j{t)) curve as an "experimental fact". This is illustrated in Fig.|Hlwhere we plot 
the curve (V^, Jmin(^)) together with projections of a selected set of AP trajectories. The AP 
solutions were obtained for a space-clamped version of 0] with initial conditions for j and V as 
shown in the figure and all other variables in their resting states. These trajectories allow us 
to follow the correlation between the transient of j and the AP V. Indeed, in the tail of an AP 
solution, the curve j vs V is almost independent of the way the AP is initiated. As a result, the 
projections of the trajectories {V{t),j{t)) intersect the critical curve {Va, jmin{Va)) in a small 
vicinity of one point, = (0.2975 ±0.0015, — 72.5±0.5). This result suggests the following 

interpretation. As a wave front propagating into the tail of a preceding wave reaches a point in 
the state corresponding to this "absolute refractoriness" point it will stop because of 

insufficient excitability of the medium, and dissipate. 

In a broader context, in the front propagation speed c is a function of j and V in the relative 
refractoriness region of the {V,j) plane, so the highly correlated dependencies of V{t) and j{t) 
in the wake of an AP mean that c at a particular point becomes a fixed function of time. 
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This makes it possible to describe c in terms of the diastohc interval (DI), i.e the time passed 
after the end of the preceding AP. This dependence, known as dispersion curve or velocity 
restitution curve, is an imp ortant tool in simplified analysis of complex regimes of excitation 



propagation 



4.3 Propagation block in two dimensions 

In two spatial dimensions, the condition of dissipation j < may happen to a piece of a 
wave front rather than the whole of it. In that case we observe a local block and a break- 
up of the excitation wave. Figure IHl shows how it happens in a two-dimensional simulation 
of the detailed model of Courtemanche et al. 0. A spiral wave was initiated by a cross-field 
protocol. This spiral wave develops instability, breaks up from time to time, and eventually 
self-terminates. This is one of the simulations discussed in detail in 0. Here we use it to test 
our newly obtained criterion of propagation block. The red colour component represents the V 
field, white for the resting state and maximum for the AP peak. This is superimposed onto an 
all-or-none representation of the j field, with black for j > and blue for j < j*. Thus the red 
rim represents the 'active front' zone where excitation has already happened but j gates are 
not de-activated yet; most of the excited region is in shades of purple representing the gradual 
decay of the AP with j deactivated. The wave ends up with a blue tail, which corresponds to V 
already close to the resting potential but j not yet recovered and still below j=K . So the blue zone 
is where there is no excitation, but propagation of excitation wave is impossible, i.e. absolutely 
refractory zone. The black zone after the tail and before the new front is therefore relative 
refractory zone, where front propagation is possible. Thus, in terms of the colour coding of 
figure IHl the prediction of the theory is: the wavefront will be blocked and dissipate where and 
when it reaches the blue zone, and only there and then. This is exactly what happens in the 
shown panels: the red front touches the blue tail, first at the third panel, at the point indicated 
by the white arrow, and subsequently in its vicinity. The excitation front stops in that vicinity 
and dissipates. So we have a break-up of the front. 

The analysis of the numerics, which ran for the total of 7400 ms until self-termination of 
the spiral and showed 4 episodes of front break-up, has confirmed that in all cases the break-up 
happened if and only if the front reached the blue region j < j*. 



4.4 Curvature effects 

Since we attempt to compare the results of our one-dimensional model to simulations of spiral 
waves in two-dimensions, it is important to explore the dependence of the solution on the 
curvature of the front. The standard theory says that in two dimensions the normal velocity of 
the wave front need to be corrected by the term A/C where A is the typical width of the wave front 
[2^. The speed-curvature diagram presented in Fig. llOf a) shows that in our simplified model 
this relationship is satisfied to rather large values of |/C|. Our choice of boundary conditions in 
Eqs. IHl assumes that the excitation fronts propagate from right to left, so positive values of the 
curvature correspond to concave fronts. Only at very small values of the radius of curvature of 
the order of 0.3 mm for j = 1 the wave speed shows a non-linear dependence on curvature as 
seen in the insert Fig. llOf b). This part of the figure also demonstrates that there is a critical 
value of the curvature for which the excitation wave stops to propagate as well as an unstable 
branch of the solution. However, these phenomena occur at very large curvatures which are far 
outside of the range of values of |/C| < O.lmm^^ observed in the two-dimensional simulations of 
Fig. El 
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The most important question with respect to om' study is whether the curvature changes 
significantly the critical value of the excitation parameter j^, below which the wave fronts fail 
to propagate. To answer this question we present Fig. llOf c) in which the wave speed c is shown 
as a function of j for three values of the curvature corresponding to a non-curved front and to 
convex and concave fronts with radius of curvature equal to 10 mm. The values of jmin for these 
three cases differ only slightly. So, the propagation blocks in our simulations do not depend 
significantly on the curvature of the front. 

This conclusion is valid for the particular cardiac model 0] and for the particular context. 
In [i^ . the minimal diastolic interval, defined as time from the moment V = — 50mV to 
the moment propagation becomes possible again, depended only slightly on curvature for the 
Modified Beeler-Reuter model at standard parameters, but was much more pronounced when 
Tj was artificially increased 6-fold. The simplest explanation of this difference is that the small 
variation of jmin due to the curvature takes much longer for j{t) to make if dj/dt is very small, 
so even that small variation jmin becomes significant. 

5 Conclusions 

In this paper, we have shown that propagation of excitation and its block in Courtemanche 
etal. a model of human atrial tissu; can be successfully predicted by a simplified model of the 
excitation front, obtained by an asymptotic description focussed on the fast sodium current, 
-^Na- Whereas it was known that main qualitative features of the iNa-driven fronts can be 
described by a two-component model for V and h, we have now found that for good quantitative 
predictions, one must also take into account the dynamics of m gates. Thus, we have proposed 
a three-component description of the propagating excitation fronts given by Eqs. ID We have 
obtained an exact analytical solution for piecewise-linear 'caricature' three-component model of 
Eqs.^ For an appropriate choice of parameters, it reproduces the key qualitative features of the 
accurate three-component model of Eqs.0]and gives a correct order of magnitude quantitatively. 
Numerical solution of the automodel equation of the proposed three-component model of Eqs.0 
gives a very accurate prediction of propagation block in two-dimensional re-entrant waves. For 
the given model, this reduces to a condition involving the pre- front values of V and j, or even in 
terms of j alone. This provides the sought-for operational definition of absolute refractoriness 
in terms of j, simple and efficient. 

The success of the propagation block prediction justifies the assumptions made on the asymp- 
totic structure, i.e. appearance of the small parameter e, of Eqs. |21 and also confirms that 
two-dimensional effects, e.g. front curvature, do not significantly affect the propagation block 
conditions, at least in the particular simulation. 

As the description and role of /Na are fairly universal in cardiac models, most of the results 
should be applicable to other models. However, some other cardiac models may require a more 
complicated description. For instance, the contemporary 'Markovian' description of /Na 
e.g. ] is very different from the classical m^jh scheme. Also, propagation in ventricular tissue 
in certain circumstances can be essentially supported by L-type calcium current rather than 
mostly /Na alone [5ll |. 

A Numerical method 

For a numerical solution the problem needs to be formulated on a finite interval z G [zmim Zmax] rather 
than on the open interval z € (— oo, oo). Furthermore, because of the piecewise definition of the problem 
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this interval must be separated in three parts [z„iin, 0], [0, £] and z,-nax] as discussed in the beginning of 
scction|21 The standard numerical methods we use require that the problem is posed on a single interval, 
for instance y G [0, L]. So we use the mapping, 

{-z, z e [z„ii„, 0], 

zG[0,e], (19) 

to transform Eqs.Has follows 

K=-{cTh{Vi)y\l-hi), 
m[ = (cT,„(yi)) ^mi, 

V^' = ((c - k) Vi - gNa{VNa ~ V2) J mf)/p, 

m'2 = -(pCTm{V2)) ^m2, (20) 
VI' = [c~k) VI - gNa{VNa " V^) J /ig m^, 

h', = -{cTH{V3)y\3, 

m'^ ^ {cTrniVs)) ^(l-TOa), 

c = 0, 

p' = 0, where p = ^/L 

V' =0 

where the subscripts 1, 2 and 3 denote the variables corresponding to the three subintervals. Here, 
the end of the second subinterval ^ is an unknown parameter and together with the wave speed c and 
the post-front voltage Kj must be determined as a part of the solution. Because these unknowns are 
constants, their derivatives must vanish which leads to the introduction of the last three equations in 
Eqs.lSOl 

The boundary conditions in Eqs.|Slat infinity are substituted by 



= [u\ +v, (21) 

' (too) 



where u is the vector of unknown variables and is a vector of small perturbations, obtained as a solution 
of Eqs. [3 linearised about Eqs.|Hl Together with the implicit assumptions V^(0) = Vm and V{^) = Vh 
which break the translational invariance and the additional requirements that the solutions must be 
continuous functions of z and that V{z) must be smooth, the necessary 15 conditions are 

^1 (0) = Vh , V2{Q)=Vh, ^3 (0) = K^ , 

Vi{0) = -p{Q) Vi{Q), hi{0) = h2{0), TOi(O) = m2(0), 
^3'(0) = p{L) V^{L), hsiO) = h2{L), TO3(0) = m2{L), 

Vi{L) = -{c{L) - {Vi{L) + 14), V2{L) = (22) 
Vm = -{V3{L) - V^L))/ {c{L)th{Vs{L))) , 
hi{L) = l, TOi(i)=0, 

^'^^^ = ■Z'^^\,(T^\ i^7V^lVrn\ + - '^^l ' 

gNaJ (yNa-V3[L)) ^ c(L) T/^ (V3 (i)) J 

We use the boundary-value problem solver D02RAF of the NAG numerical library which emplo ys a finite- 
difference discretization coupled to a deferred correction technique and Newton iteration |53|. The 
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analytical solution given in Eas. 1101 is used as an initial approximation to start the correction process. 
The method proves to be very robust over a large range of parameters. 

The authors are grateful to I.V. Biktasheva for sharing her experience of simulation of model 0|, 
to H. Zhang and P. Hunter for inspiring discussions related to this manuscript and to the anonimous 
referees for constructive criticism and helpful suggestions. This work is supported by EPSRC grants 
GR/S43498/01 and GR/S75314/01. 
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Model 


Wave speed 


Rel. error 


Post-front 
voltage 


Maximum rate 
of AP rise 




C, [mm/ms] 


in C 


K., [mV] 


(dy/dt)„,ax [v/s] 


The full model of 










Courtemanche et al. [6] 


0.2824 




3.60 


173.83 


Model [6] with replacements 










h{v)^h{v)e{Vh-v), 










m{V) -^fn{V)e{V -Vm) 


0.2130 


24.5 % 


-0.99 


173.83 


Equations |31 with 










M{V) =m{V), H{V) = h{V) 


0.2095 


25.8 % 


-1.06 


183.82 


Equations |21 with 










M{V) = 1, H{V) = 1, i.e. Eqs.H 


0.2372 


16.0 % 


2.89 


193.66 


Equations El 


0.4422 


57.3 % 


18.26 


643.97 



Table 1: A comparison of the wave speed C and post-front voltage amplitudes V^j and the 
maximum rate of AP rise (dV/ dt)max of various approximations to the Courtemanche et al. 
0] model. Prior to firing, the tissue in the models was set at rest at the standard values of the 
parameters, see 0|. In these and other numerical results /C = is assumed unless explicitly 
stated otherwise. Space-clamped versions of the models are used to compute ( dV/ di)max- 
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Figure 1: Asymptotic properties of the atrial model of Courtemanche et al. (a) Time scale 
functions of dynamical variables vs. time, (b) Quasistationary values of the gating variables m 
and h. (c) Transmembrane voltage as a function of time, (d) Main ionic currents vs. time. 

-^in = -^b,Na + -^NaK + -^Ca,L + -^b,Ca + -^NaCa and /out = /p,Ca + -^Kl + ^to + ^Kur + -^Kr + ^Ks + ^b.K 

are the sums of all inward and outward currents, respectively and the individual currents are 
described in 0]. The results are obtained for a space-clamped version of the model at values 
of the parameters as given in |^. In (c) and (d) a typical AP is triggered by initialising the 
transmembrane voltage to a non-equilibrium value V = —20 mV. 
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Figure 2: Response to a temporary local block of excitability (B) in the models of (A) Courte- 
manche et al. ja], (B) FitzHugh-Nagumo Eqs.nand (C) in Eqs. 01 The border of the blocked 
region is shown by broken lines. Solutions are represented by shades of gray: black is the smallest 
and white is the largest value of V within the solution. The parameters of the FitzHugh-Nagumo 
model are f3 = 0.75, 7 = 0.5 and €g = 0.03, while for the two other models the same parameter 
values as described in 0| are used, the block is described in the plots. The value of j = 0.28 in 
the block in (c) is just below the propagation threshold, see Fig. |H1 The time and space ranges 
(in dimensionless units) are 70 x 70 in (B) and 80 x 50 in (A) and (C). 
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Figure 3: (A) The AP potential and (B) the gating variables h and m as functions of the 
travelling wave coordinate Z = z\fD. The solution of the model of Courtemanche et al. 0] is 
given by circles, of the full three- variable model of Eqs.^by thin lines, and the analytical solution 
given by Eqs.UDlfor TnI = l^^^Vm) = 781.8, = ThiVm) = 1.077, = r^(Ki) = 0.131, 
Va = —81.18 mV and j = 0.956 by thick lines. The gates h and m are indicated in the plot. 
The position of the internal boundary point H = £,^fD is indicated by a dash-dotted line. 



Conditions for propagation in atrial tissue 



23 



0.4 
^0.3 



0.1 



" 0.2 0.4 0.6 0.8 1.0 

j 

Figure 4: The wave speed C as a function of the excitation parameter j. Thick Hnes: the 
numerical solution of Eqs. [7| Thin lines: solution Eq. E]for values of t/j and Tm. corresponding 
to a selected voltage V = in Eqs. El From right to left: K = —28, —30, Vm, —34, —36 
—38 (mV). In both cases Va = —81.18 mV and /C = mm~^. 
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Figure 5: The wave speed C as a function of the time-scale ratio Th/Tm in the caricature 
model Eas. I7I8I The values of Th and /^a are fixed to the values of the corresponding functions 
in Eqs-ISlat a selected voltage F = K, the pre-front voltage is Va = —81.18 mV and curvature 
is /C = mm-^ Left plot: left to right, K = -38, -36, -34 and -32.7 = Vm (mV), and 
j = 0.9775. Right plot: right to left, j = 0.2 to 1.0 and K = Vm- 
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Figure 6: The threshold value jmin above which propagation is possible, as a function of the pre- 
front voltage Va for the same values of the parameters as in Fig.|31 i.e. Th = 1.077, = 0.131. 
Shown are different approximations to the perturbation expansion given by Eq. 1161 solid line, 
zeroth order, Eq. 1171 dashed line, first-order, Eq. 1181 dotted line: second-order. 
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Figure 7: The wave speed C as a function of j and Va, for the model of Eqs.[7| Rapid changes 
are indicated by a higher density of curves. The thick dotted hne on the base represents the 
threshold value jmin and may be compared to the results in Fig. El 
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Figure 8: The thick sohd hne represents the threshold value jmin for excitation failure as a 
function of Va for the model given by Eqs. [3 The dotted lines represent projections of AP 
trajectories in the space-clamped detailed model of 0. 
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Figure 9: Local propagation block, dissipation and break-up of the front of a re-entrant exci- 
tation wave. The density plots represent the distribution of the transmembrane voltage V (red 
component) in regions of super-threshold (white) and of sub-threshold (blue) excitability j. The 
white arrow indicates the time and place the propagation block begins. The time increases from 
(A) to (F) with At = 20 ms; size of the simulation domain is 75 mm x 75 mm. 
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Figure 10: (A) and (B) The wave speed C for the model of Eqs. [3 |S1 as a function of the 
curvature for values of j = 1 . . . 0.4 (from top to bottom). Results for the detailed model 0| are 
denoted by thick solid lines. (C) The wave speed C in the model given by Eqs.[7|as a function 
of j for /C = 0.1, and —0.1 mm"'^ (from top to bottom). 



